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We present an improved variant of the in-medium similarity renormalization group (IM-SRG) 
based on the Magnus expansion. In the new formulation, one solves flow equations for the anti- 
hermitian operator that, upon exponentiation, yields the unitary transformation of the IM-SRG. 

The resulting flow equations can be solved using a hrst-order Euler method without any loss of 
accuracy, resulting in substantial memory savings and modest computational speedups. Since one 
obtains the unitary transformation directly, the transformation of additional operators beyond the 
Hamiltonian can be accomplished with little additional cost, in sharp contrast to the standard 
formulation of the IM-SRG. Ground state calculations of the homogeneous electron gas (HEG) and 
nucleus are used as test beds to illustrate the efficacy of the Magnus expansion. 
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I. INTRODUCTION 

The quest to predict and understand the properties 
of exotic nuclei starting from the underlying nuclear 
forces represents a cornerstone of modern nuclear the¬ 
ory. Already for stable nuclei, there are computational 
and theoretical challenges that make the ab-initio de¬ 
scription of nuclear structure quite difficult. Neverthe¬ 
less, tremendous progress has been made over the past 
two decades, where it is now possible to perform quasi- 
exact calculations including three-nucleon interactions of 
nuclei up through Carbon or so in Quantum Monte Carlo 
(QMC) and No-Core Shell Model (NCSM) calculations, 
and N = Z nuclei up through ^®Si in lattice effective field 
theory with Euclidean time projection [IH3]. 

Since exact methods scale unfavorably with system 
size, it is necessary to develop approximate, but sys¬ 
tematically improvable methods to extend the reach of 
ab-initio theory beyond light nuclei. Over the past 
decade, Coupled Cluster (CC) theory, Self-Consistent 
Green’s Functions (SCGF), Auxiliary Field Diffusion 
Monte Garlo (AFDMG), and the In-Medium Similar¬ 
ity Renormalization Group (IM-SRG) have been success¬ 
fully applied to calculate properties of selected medium 
mass nuclei and infinite nuclear matter [iHin]- Early 
applications of these methods were limited primarily to 
ground state properties of stable nuclei near shell closures 
with two-nucleon forces only. In recent years, however, 
substantial progress has been made on including three- 
nucleon forces [giiiiiiii], targeting excited states and 
observables besides energy [laiii], and moving into the 
more challenging terrain of open-shell and unstable nu¬ 
clei [TsHT^ . 

The IM-SRG is a particularly appealing method due 
to its flexibility to target ground and excited state prop- 
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erties for closed- and open-shell systems. As discussed in 
Sec.jnj the essence of the IM-SRG is to perform a contin¬ 
uous unitary transformation on the Hamiltonian (and all 
other observables of interest) to drive it to a diagonal or 
block-diagonal form. The transformation is implemented 
by solving a coupled set of flow equations for the matrix 
elements of the Hamiltonian and any other operators of 
interest 

H{s) = U\s)HU{s) ^ = [vis),H{s)] 

0{s) = UHs)OU{s) ^ ^ = [r;(s), 0 (s)], 

( 1 ) 

where s is a continuous flow parameter, and the choice of 
the generator r]{s) = implicitly defines the trans¬ 

formation U{s). Despite the flexibility to tailor 77 to a 
wide range of problems and the modest computational 
scaling with system size, the formulation in Eq. suffers 
from the following difficulties: 

• The coupled ODEs can become stiff for certain 
choices of generator and/or for systems with strong 
correlations. 

• The numerical integration of Eq. requires a high- 
order ODE solver to accurately preserve the eigen¬ 
values of the evolved Hamiltonian. The use of a 
high-order solver consumes a large amount of mem¬ 
ory since multiple copies of the solution vector (e.g., 
15-20 for the predictor-corrector solver of Shampine 
and Gordon m) need to be stored at each time 
step. 

• For each additional observable of interest, the num¬ 
ber of coupled ODEs that need to be solved is 
roughly doubled, assuming a comparable level of 
truncation for the evolved operator as the Hamil¬ 
tonian. Moreover, the flow equations for the ad¬ 
ditional observable(s) can exacerbate the problems 
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with stiff ODEs, since the time scales for the oper¬ 
ator evolution may be very different from those of 
the Hamiltonian. 

In the present paper, we will demonstrate how these 
difficulties can be circumventec03y using the Magnus ex¬ 
pansion to recast Eq. as a flow equation for the opera¬ 
tor H(s), where U{s) = The unitary operator C/(s) 

is subsequently used to transform the Hamiltonian and 
any other operators of interest via the Baker-Cambell- 
Hausdorff (BCH) formula. We will show that in the Mag¬ 
nus expansion formulation, one can use a naive hrst-order 
forward Euler method to solve the flow equations for f2(s) 
without any loss of accuracy. This provides a substantial 
reduction in memory consumption, and allows operators 
beside the Hamiltonian to be evolved with little addi¬ 
tional cost, in sharp contrast to the original formulation 
of the IM-SRG based on direct integration of Eq. 

The remainder of the paper is organized as follows: 
In Section |llj we review the basic formalism of the SRG 
and illustrate how the Magnus expansion can be used 
to make its numerical implementation more efficient for 
a schematic model. In Section |HI[ we give some imple¬ 
mentation details of our IM-SRG and Magnus expansion 
calculations of the homogeneous electron gas (H EG) and 
nucleus. Results are presented in Section |IV[ and 
conclusions are presented in Section [V| 

II. FORMALISM 
A. SRG 

The similarity renormalization group consists of a con¬ 
tinuous sequence of unitary transformations that grad¬ 
ually suppress off-diagonal matrix elements, driving the 
Hamiltonian towards a band- or block-diagonal form [ 221 - 
[25] . Writing the transformed Hamiltonian as 

H{s) = U{s)HU^ (s) = H^{s) + (2) 

where H‘^{s) and are the arbitrarily defined “di¬ 

agonal” and “off-diagonal” parts of the Hamiltonian, the 
evolution with the continuous flow parameter s is given 
by 

^ = [,7(s),7I(s)], (3) 

where the 77 ( 5 ) = U{s)dW (s)/ds is the (anti-hermitian) 
generator of the transformation. The choice of the gen¬ 
erator first suggested by Wegner, 


^ See also the recent Driven Similarity Renormalization Group 
(DSRG) approach of Ref. |21| . where the problem is recast so 
that one solves non-linear amplitude equations instead of flow 
equations. 


guarantees 

^Tr = 2Tr{rf) = -2Tr{r)'<r]) < 0, (5) 

ds 

which demonstrates that the strength of decays 
with increasing s |23j . By analyzing the flow equations 
in the eigenbasis of H‘^{s) and defining Hf^{s) = e^, 
one can show that the Wegner generator gives a super¬ 
exponential decay of the off-diagonal matrix elements 

^ . ( 6 ) 

The SRG evolution with the Wegner generator closely 
resembles the conventional Wilsonian RG, since ma¬ 
trix elements between widely-separated energy scales are 
eliminated before moving inwards towards the diagonal. 
The Wegner generator is numerically very stable, but the 
different rates of decay for off-diagonal matrix elements 
can lead to stiff ODEs. To avoid this, two alternative 
classes of generators are commonly used in nuclear ap¬ 
plications. The first alternative was proposed by White 
in Ref. [26] 



which leads to a uniform suppression of off-diagonal ma¬ 
trix elements 

Ht^{s) ^ e-^H°f{Q). ( 8 ) 

The White generator is numerically very efficient for well- 
behaved problems, though it can become unstable when 
the energy denominator in Eq. becomes too small. In 
Ref. [7|, the so-called “imaginary time” generator was 
proposed as a compromise between the White and Weg¬ 
ner generators 

77„(s) = sign{€i - ej)H°f{s ), (9) 

where the leading behavior of the off-diagonal matrix el¬ 
ements is 

~ . (10) 

In the present paper, all of our IM-SRG calculations were 
done using the White generator. However, we stress 
that the computational benefits of the Magnus expansion 
carry over irrespective of the specific choice of generator. 

B. In-medium Evolution 

Until recently, most SRG applications to nuclear in¬ 
teractions have been carried out in free space to “soften” 
two- and three-nucleon interactions to be used as input 
for ab-initio calculations EilllZlISH]. The free-space evo¬ 
lution is convenient, as it does not have to be performed 
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for each different nucleus or nuclear matter density. How¬ 
ever, it is necessary to consistently evolve three-nucleon 
(and possibly higher) interactions to be able to soften 
the interactions significantly and maintain approximate 
s-independence of H ^ 3 observables. The consistent 
SRG evolution of three-nucleon operators represents a 
significant technical challenge that has only recently been 
solved in recent years P5H5T] . 

An interesting alternative is to perform the SRG evo¬ 
lution in-medium (IM-SRG) for each A-body system of 
interest [SlEH]. Unlike the free-space evolution, the IM- 
SRG has the appealing feature that one can approx¬ 
imately evolve 3, ...,A-body operators using only two- 
body machinery by normal-ordering with respect to an 
A-body reference state. Moreover, with a suitable defi¬ 
nition of the off-diagonal part of the Hamiltonian to be 
driven to zero, the IM-SRG can be used as an ab-initio 
method in and of itself, rather than simply to soften the 
Hamiltonian as in the free-space SRG. 

Starting from a general second-quantized Hamiltonian 
with two- and three-body interactions. 


H — 'y ( TgrO^jUr + 2j2 ^ ^ 


qrst 


^ X! ^qfituva\alalayauat 


( 11 ) 


qrstuv 


all operators can be normal-ordered with respect to a 
finite-density Fermi vacuum |$) (e.g., the Hartree-Fock 
ground state), as opposed to the zero-particle vacuun 0 
Wick’s theorem can then be used to exactly write H as 


H — E y ) fqj- . • F . y ) ^qrst ■ a^a^Ufas • 


qr 


qrst 


H” ^ ^ qrstuv • . , 


( 12 ) 


qrstuv 


where strings of normal-ordered operators obey the fol¬ 
lowing relation. 


(<l>|:al...a, :|<1>)=0, (13) 


^ In the present work, we restrict our attention to single reference 
(i.e., closed-shell) systems for which a single Slater determinant 
provides a reasonable starting point. See Refs. [mnzuis] for 
extensions of the IM-SRG to open-shell systems. 


and the terms in Eq. (12 1 are given by 

E = TqqUq + - 


qr 


1 


F g y ) j 

qrs 

1 ^ ’ r( 3 ) 


fqr — Tqr F '^^Vqsrs''^s F ^ ^qstrst'^a'<^t > 

s st 

y-y(3) 


r — 

i qrst — ''qrst 


w _ 1/^^) 

I'y qrstuv — y qrstuv ■ 


(14) 

(15) 

(16) 
(17) 


Here, the initial n-body interactions are denoted by 
and Uq = 0(eF — Cq) are occupation numbers in the refer¬ 
ence state 1$), with Fermi energy ep. It is evident that 
the normal-ordered 0 -, 1 -, and 2 -body terms include con¬ 
tributions from the three-body interaction through 
sums over the occupied single-particle states in the refer¬ 
ence state 1$). Neglecting the residual three-body inter¬ 
action leads to the normal-ordered two-body approxima¬ 
tion (N02B), which has been shown to be an excellent 
approximation in many nuclear systems pun Eg. Trun¬ 
cating the in-medium SRG equations to normal-ordered 
two-body operators, which we denote by IM-SRG(2), will 
approximately evolve induced three- and higher-body in¬ 
teractions through the nucleus-dependent 0-, I-, and 2- 
body terms. 

Using Wick’s theorem to evaluate Eq. Pwith H{s) = 
Eo{s) + f (s) + r{s) and 77 ( 5 ) = (s) + rj^ {s) truncated 

to normal-ordered two-body operators, one obtains the 
coupled IM-SRG (2) flow equations 


— y ) Vqrfrqi’l^q ^r) F ^ ^ ^^ 'Cjqrst^ stqr'^^q'^r^^s'^t ■ 


qrst 


(18) 


df, 


ds Pqr')^qs f sr 

s 

^ ^ (^5 ’^t'jiVst^tqsr fstVtqsr) 

st 

+ '^{nsritflu + nsntnu){l + Pqr)r]uqst^ stur 


(19) 


dr, 


qrst 


ds = E{(i- Pqr^{j]qu^ urst fqu^ursty\ 

u 

^ ^ {(1 Pst^i'Hus^qrut fusVqrut)} 

u 

^ ^v^i^qruv^uvst ^qruvVuvst') 


^ ^ Pst^'^vrut^ 


uqvs ? 


( 20 ) 
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FIG. 1. Schematic representation of the initial and final 
Hamiltonians, H{0) and H{oo), in the many-body Hilbert 
space spanned by particle-hole excitations of the reference 
state. 


where fir = 1 — Ur and the s-dependence has been sup¬ 
pressed for brevity. 

For the calculation of the ground state of a closed-shell 
system in the IM-SRG(2) approximation, it is simple to 
identify fait+h.c}, where a,b denote par¬ 

ticle (unoccupied) and i,j hole (occupied) single particle 
states, as the relevant vertices which connect our chosen 
reference state |$) with higher particle-hole excitations, 
see Fig. By designing a generator to eliminate these 
terms, one finds that the 0-body term approaches the 
interacting ground state energy in the limit of large s, 

\im Eo{s) = {m{sm = Eg,. (21) 

s-^-oo 


In the present paper we use the White generator [26] . 
though we deviate slightly from recent implementations 
that use Epstein-Nesbett energy denominators [12], opt¬ 
ing for the simpler Mpller-Plosset energy denominators 


= E 


fa 


fa - fi 


E 

abij 


■ oloi 


- abij 


fa + fb- fi- fj 


- : alalajai : —H.c., (22) 


where fa = faa, etc. The use of Mpller-Plosset denomi¬ 
nators has minimal impact on the results of ground state 
IM-SRG(2) calculations, but it has the virtue of reveal¬ 
ing the connection to MBPT. In a subsequent work, this 
connection will be used to develop approximations that 
go beyond the IM-SRG(2). 


C. Magnus Expansion 


State-of-the-art solvers can require the storage of 15-20 
copies of the solution vector in memory, which becomes 
problematic for large model spaces. The problem is ex¬ 
acerbated if one wants to calculate additional observ¬ 
ables, roughly doubling the memory requirements assum¬ 
ing the same N02B truncation as for the Hamiltonian. 
Moreover, the additional flow equations for each observ¬ 
able can evolve with rather different timescales than the 
Hamiltonian, which increases the likelihood of the ODEs 
becoming stiff. 

To bypass these limitations, we now describe an al¬ 
ternative method to solving Eq. (§ using the Magnus 
expansion [32] . In the notation of our present problem, 
our starting point is the differential equation obeyed by 
the unitary transformation, 

^=-r,is)Uis), (23) 

where 17(0) = I and W{s)U{s) = U{s)W{s) = 1. 
This can be formally integrated and written as the time- 
ordered exponential 

U{s) =Ts{e-fo^^^'^‘^^'} (24) 

= 1 - ds'r]{s') + r ds' r ds''r,{s')r]{s") + ... 

Jo Jo Jo 

(25) 

Eq. is not very useful in practical calculations since 
i) there is no guidance on how the series should be trun¬ 
cated, ii) one would need to store rj for multiple s-values, 
and iii) it is not obvious how to consistently transform 
the Hamiltonian and other observables in a fully linked, 
size-extensive manner with the truncated series. 

The Magnus expansion is essentially a statement that, 
given a few technical requirements on r]{s), a solution of 
the form 


U{s) = (26) 

exists, where Hf(s) = —H(s) and fl(0) = 0. In most pre¬ 
vious applications of the Magnus expansion, one typically 
expands r 2 (s) in powers of 77 ( 5 ) as 

00 

H ^ . (27) 

n—1 

Eor issues of convergence and mathematical details, see 
Refs. (mm). Gombining this with the formally exact 
derivative 


IM-SRG calculations typically use ODE solvers based 
on high-order Runge-Kutta or predictor-corrector meth¬ 
ods to solve Eq. The use of a high-order method is 
essential as the accumulation of time-step errors will de¬ 
stroy the unitary equivalence between H{s) and 77(0), 
even if no truncations are made in the flow equations. 



fc=0 


ad^iv) = V 
adai-n) = [n,ad'^-\rj)] , 


( 28 ) 
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where Bk are the Bernoulli numbers and ad^{r]) the re¬ 
cursively defined nested commutators, one can obtain ex¬ 
plicit expressions for the 0 „(s), 

ni(s) = -/ dsi-q{si) 

Jo 

^^ 2 ( 5 ) = ^ [ dsi [ ds2[77(si),ry(s2)] (29) 

^ Jo Jo 


As expected, rewriting the time-ordered exponential as a 
true matrix exponential moves the complications of time 
ordering into the expression for r2(s). The utility of the 
Magnus expansion lies in the fact that, even if is trun¬ 
cated to low-orders in ry, the resulting transformation in 
Eq. 26 using the approximate is unitary, in contrast to 
any truncated version of Eq. |24[ 

For large-scale IM-SRG calculations, the expressions in 
Eq. are of limited value since they require the storage 
of 77 ( 5 ) over a range of s-values. Therefore, in the present 
work we instead construct r 2 (s) by numerically integrat¬ 
ing Eq. subject to certain approximations discussed 
below. The transformed Hamiltonian, and any other op¬ 
erator of interest, can then be constructed by applying 
the Baker-Cambell-Hausdorff (BCH) formula, 


H{s) = e^He-^ = j:Zo ( 30 ) 

O(s)=e^Oe-^=Er=0M“4(O). (31) 

Before discussing how we truncate Eqs. and in 
practical calculations, it is instructive to study a simple 
matrix model that can be solved without any truncations. 
Consider the initial Hamiltonian 

i? = T + E= , (32) 

where the diagonal “kinetic energy” term is 

T = (J f 1 ) . (33) 


Let us now try to diagonalize H using the Wegner gen¬ 
erator T]{s) = [T, 7J(s)], solving the SRG equations using 
the Magnus expansion and by direct integration of Eq. 
Note that for this choice of initial H, both r]{s) and r2(s) 
are real, antisymmetric matrices throughout the flow 


vis) = i9rjis)(J2 (34) 

V,{s) = ign{s)a 2 , (35) 


where a 2 is the Pauli matrix. Consequently, Eq. [^termi¬ 
nates at the first term and Eq.j^can be summed up to all 
orders using the well-known properties of Pauli matrices. 
Since the large memory footprint of high-order adaptive 
solvers is the main computational challenge in large-scale 
SRG calculations, let us instead try to use a naive first- 
order Euler method to integrate Eqs. [^ and 28 The 



s 


FIG. 2. (Color online) |4/ii(s) — Sgsl versus s for different 
Euler step sizes calculated via direct integration of the SRG 
flow equation, Eq.j^ and using the Magnus expansion, Eqs.[28| 
and[30l 


results are shown in Fig. [^ where we plot |iLii(s) — Egs\ 
- which should go to zero at large s — versus s for differ¬ 
ent Euler step sizes Ss. Unsurprisingly, we see that the 
direct integration of Eq. accumulates large time-step 
errors, with the plateaus at large s displaying a strong 
dependence on the Euler step size. The Magnus solution, 
on the other hand, converges to a final answer at large s 
that is independent of step size and agrees with the exact 
result to within machine precision. The insensitivity to 
the time ste p si ze is due to the fact that while each Euler 
step in Eq. 28 gives an error of order C)((5s^), the ex¬ 


ponentiated operator at the end of the evolution is still 
unitary. This is the primary advantage of the Magnus 
expansion; by reformulating the problem to solve flow 
equations for H(s) instead of H{s), one can use a simple 
first-order Euler method and dramatically reduce mem¬ 
ory usage. Once r2(s) is in hand, the transformation of 
His) and any other observables of interest immediately 
follows from Eq. In contrast to the direct integration 
of Eq. the dimensionality of the flow equations does 
not increase when one evolves additional observables. 


D. Magnus(2) Approximation 

Having illustrated the advantages of the Magnus ex¬ 
pansion in a simple model, we would now like to ap¬ 
ply it in many-body calculations. Unlike in coupled 
cluster theory where the BCH formula for the similar¬ 
ity transformed Hamiltonian terminates at finite order, 
Eqs. [^ and [^ involve an infinite-order series of nested 
commutators that generate up to A-body operators. To 
make progress, we introduce the Magnus(2) truncation in 
which all commutators (as well as r 2 (s), 77 (s) and His)) 
are truncated to the N02B level. Even with this ap¬ 
proximation, the expressions for dVL/ds and iL(s) involve 
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an infinite number of terms. However, for both Eqs. 28 
and [30| at the N02B level, we empirically find that the 
magnitude of terms decreases monotonically in k for all 
systems studied thus far. Therefore, we numerically trun- 
at the term if 


cate Eqs. 


28 


Bk\\ad^{ri)\\ 


fc!||H|| 


^deriv ■ 


(36) 


For the truncation of ( |30[ ) , we could use a similar crite¬ 
ria as for the derivative expression. However, since we are 
interested in the ground-state energy, we use a simpler 
condition where the series is truncated when the zero- 
body piece of the k*"^ term falls below some threshold. 


WdUH)}ob 

k\ 


< eBCH ■ 


(37) 


In the calculations presented below, we will find that the 
final results are insensitive to large variations in Cderiv 
and ebch, which we take as an a posteriori justification 
for our truncations. 


III. HAMILTONIANS AND IMPLEMENTATION 


Before presenting the results of IM-SRG(2) and Mag- 
nus(2) calculations of the homogeneous electron gas 
(HEG) and ^®0, we review some details of our imple¬ 
mentations for both systems. For the homogeneous elec¬ 
tron gas, we perform our calculations for the closed-shell 
configuration of iV = 14 electrons in a cubic box with pe¬ 
riodic boundary conditions. Note that if one is interested 
in extrapolating to the thermodynamic limit, calculations 
should be done for a larger closed-shell configurations of 
N = 38, 54, 66,... electrons, with finite-size corrections 
for the kinetic and potential energy taken into account. 
Here we neglect these corrections since our primary pur¬ 
pose is to demonstrate the effectiveness of the Magnus 
expansion, and the quasi-exact Full Configuration Inter¬ 
action Quantum Monte Carlo (FCIQMC) results we com¬ 
pare against also neglect these corrections [3S]. The rel¬ 
evant single particle orbitals are plane waves with quan¬ 
tized momenta 

^k.(r) = ^e*‘^'-X., (38) 


where is the box volume, x<y is a spin eigenfunction, 
and k = ^(n^;, Uj,, n^) where rix^riy, and are integers. 
We follow common practice and use the Wigner-Seitz ra¬ 
dius to characterize the density of the HEG, 


rs = 


ao ’ 


(39) 


where ag is the Bohr radius and rg is defined in terms of 
the inverse density as 



N ■ 


(40) 


We use a basis set truncation which keeps M single parti¬ 
cle states with energy less than some cutoff Ec, although 
other choices are possible [5B] . 

In the plane wave basis, the kinetic energy matrix ele¬ 
ments are diagonal 

and the Coulomb matrix elements are given by 

^ijkl — kfc<^q,ki—kj • (42) 

Note that the q = 0 term is omitted due to its can¬ 
cellation against the inert, uniform positively charged 
background that is needed to make the system charge 
neutral m- Since we are interested primarily in the cor¬ 
relation energy, we have omitted the Madelung term in 
all of our calculations. 

For the calculations of ^®0, our starting point is the 
intrinsic nuclear A-body Hamiltonian 

77 = ^ ( 43 ) 

where is the two-body part of the intrinsic kinetic 
energy, and we restrict our attention to two-nucleon in¬ 
teractions only. Results are presented for input NN in¬ 
teractions derived from the N^LO (500 MeV) potential of 
Entem and Machleidt [38] at several different free-space 
SRG resolution scales, A = 2.0, 2.7, and 3.0 fm“^. 

For both systems, the Magnus(2) and IM-SRG(2) cal¬ 
culations start by normal ordering the Hamiltonian with 
respect to the HF ground state. In the case of the HEG, 
translational invariance implies the HF orbitals are plane 
waves. Therefore, the HF reference state is just a Slater 
determinant comprised of the lowest energy doubly oc¬ 
cupied plane wave states for TV = 14 electrons. For 
we must self-consistently solve the Hartree-Fock equa¬ 
tions by expanding the unknown HF orbitals in a har¬ 
monic oscillator basis truncated to oscillator states obey¬ 
ing 2n -|- Z < Cniax; where Cmax is sufficiently large so 
that the results are approximately independent of the 
huj value of the underlying oscillator basis. For the NN 
interactions used in the present calculations, a cutoff of 
Gmax = 8 is Sufficiently large for most purposes. Once a 
converged HF ground-state is obtained, the Hamiltonian 
is normal-ordered w.r.t. to this solution, and the result¬ 
ing in-medium zero-, one-, and two-body operators serve 
as the initial values for the Magnus(2) and IM-SRG(2) 
flow equations. These are subsequently integrated until 
sufficient decoupling is achieved, as determined by the 
size of the second-order many-body perturbation theory 
MBPT(2) contribution of the flowing Hamiltonian H{s) 
to the ground state energy. We use a threshold of I0“® 
Hartree (MeV) for the HEG (^®0) calculations, respec¬ 
tively, which corresponds to relative changes in the flow¬ 
ing ground-state energy of 10“^ or less for both systems. 











7 



s s 


36 


FIG. 3. (Color online) Relative importance of the term in the Magnus derivative as defined by the lefthand side of Eq. 
evaluated in the N02B approximation. The top row is for the homogeneous electron gas at Wigner-Seitz radii of a) = 0.5 
and b) Vs = 5.0. The bottom row is for ^®0, starting from the chiral NN potential of Entem and Machleidt |38| . softened by 
a free-space SRG evolution to (c) A = 2.0 fm“^ and (d) A = 3.0 fm“^. The electron gas calculations were done for = 14 
electrons in a periodic box with M = 114 single particle orbitals. The calculations were done in an Cmax = 8 model space, 
with hui = 24 MeV for the underlying harmonic oscillator basis. 


IV. RESULTS 


We begin by examining the numerical evidence for 
truncating Eqs. and by hand. In Figure we 
plot the lefthand side of Eq. 36 for the HEG (top row) 
and (bottom row) as a function of the flow param¬ 
eter. To assess the role of correlations, the HEG cal¬ 
culations were performed at two different Wigner-Seitz 
radii, = 0.5 and = 5.0, and the calculations 
were done using NN interactions at two different resolu¬ 
tion scales, A = 2.0 fm“^ and A = 3.0 fm“^. For the 
HEG, the Ts = 0.5 contributions are completely negligi¬ 
ble by the k = 2 term, which is not surprising since the 
kinetic energy dominates in this weakly correlated high- 


density regime m- Even for the rg = 5.0 case, where 
correlations and non-perturbative effects start to become 
sizable, one finds that the successive terms in Eq. de¬ 
crease monotonically, though the individual terms are 
substantially larger than for the rg = 0.5 case. Anal¬ 
ogous results are found for the individual terms are 
larger for the harder A = 3.0 fm“^ interaction since the 
system is more strongly correlated, but they systemati¬ 
cally decrease with increasing order k. 

Figure tells a similar story for the BCH formula, 
where the lefthand side of Eq. [^is plotted as a function 
of the flow parameter. In all cases, we see the impor¬ 
tance of successive terms decreases monotonically. Reas¬ 
suringly, we find that the final results in our calculations 























are essentially independent of the convergence criteria 
provided Cderiv 10“^ and cbch ^ 10 “^, where the lat¬ 
ter is in units of Hartree (MeV) for the HEG (^®0) cal¬ 
culations, respectively. For instance, raising both con¬ 
vergence criteria from 10“® to 10“'* changes the ground 
state energy at the 1 eV (10“^ Hartree) level in the 
(HEG) calculations, respectively. 

As was illustrated for the toy model in Section |H G[ 
the key advantage of the Magnus expansion is that one 
can use a first-order Euler method to accurately solve 
the flow equations. We now demonstrate that the same 
conclusion holds for realistic IM-SRG calculations. Re¬ 
ferring to Figs. and we show the 0-body part of the 
flowing Hamiltonian H(s) versus the flow parameter for 
the electron ga^ and *®0. The black solid lines denote 
the results of a standard IM-SRG(2) calculation using 
the adaptive ODE solver of Shampine and Gordon, while 
the other curves denote IM-SRG(2) and Magnus(2) cal¬ 
culations using a first-order Euler method with different 
step sizes Ss. For the electron gas, the exact FCIQMC 
results [35] are shown for reference. Unsurprisingly, the 
IM-SRG(2) Euler calculations are very poor, with the 
various step sizes converging to different large-s limits. 
The Magnus(2) calculations, on the other hand, converge 
to the same large-s limit in excellent agreement with the 
standard IM-SRG(2) results. The insensitivity to step 
size is due to the fact that the time step errors accumu¬ 
late in D(s) as opposed to H{s). At the end of the flow, 
D(s) is still an anti-hermitian operator, and the transfor¬ 
mation in Eq. is unitary, up to truncation errors in 
the N02B approximation. 

Given that the Magnus(2) results are independent of 
step size over the range considered, one might try to keep 
increasing the step size to reach the ground state in fewer 
steps. This unfortunately is not possible, as the flow 
tends to diverge with too large of a time step. One of the 
strengths of the SRG approach is that the transformation 
is adapted as the Hamiltonian is transformed. With too 
large of a time step, we rob the method of this flexibil¬ 
ity and run the risk of applying a “large rotation” of the 
Hamiltonian that induces large three- and higher-body 
components. This would not be a problem if we evalu¬ 
ated the BCH and Magnus derivative expressions with¬ 
out approximation; the method would find its way back 
since the large rotation is still unitary if no truncations 
are made. However, in the Magnus(2) approximation we 
make, the neglect of the induced three- and higher-body 
terms can lead to a lack of convergence. Empirically, we 
find that this difficulty is avoided by enforcing that at 
each time step the “off-diagonal” matrix norm ||i?°‘*|| is 
decreasing. This can be implemented by using a simple 
mid-point integrator algorithm and decreasing the time 
step if ||i?°‘*|| has increased between the first and second 
half of a step. 


® For the HEG, we plot Eo{s) — Ehf, which approaches the cor¬ 
relation energy at large s. 


As a final demonstration of the utility of the Magnus 
expansion, we turn to the evolution of operators other 
than the Hamiltonian. In the conventional approach 
based on the direct integration of Eq. the dimension¬ 
ality of the flow equations increases with each additional 
operator to be evolved. In contrast, in the Magnus ex¬ 
pansion the dimensionality of the flow equations does not 
change; the additional computational expense shows up 
only in the evaluation of the BCH formula for the trans¬ 
formed operator, Eq. |31[ For a given operator O, we 
have 

(«'o|0|^'o) = lim , (44) 

S—>-00 

where |4)) is the reference state. Therefore, the 0-body 
piece of the transformed operator approaches the inter¬ 
acting ground state expectation value in the large-s limit. 

As a proof-of-principle, we perform a Magnus(2) evo¬ 
lution to evaluate the ground state expectation value of 
the momentum distribution operator hj. = for the 
HEG, and the generalized center of mass (COM) Hamil¬ 
tonian for the nucleus, 

i?cm(w) = Tcm -f . (45) 

Figure [7] shows the Magnus(2) ground state momentum 
distribution for a system of A^ = 14 electrons in a periodic 
box for several different Wigner-Seitz radii. Even with 
the neglect of finite size corrections and the extremely 
coarse momentum grid due to the small box sizes consid¬ 
ered, the qualitative behavior agrees with expectations 
for the electron gas; correlations become more important 
at larger rg , leading to a stronger depletion of modes with 
k < kp and smaller discontinuity at the Fermi surface. 
We note that the Magnus(2) results are in good agree¬ 
ment with the IM-SRG(2) calculations based on Eq. [^as 
well as results generated by the Feynman-Hellman theo¬ 
rem, but at a fraction of the cost. In addition to provid¬ 
ing a memory-efficient means for evolving operators be¬ 
yond the Hamiltonian, Fig. ^ shows that the Magnus(2) 
approximation gives a small but robust computational 
speedup for a range of basis sets, even including the ad¬ 
ditional effort of generating the momentum distributions, 
which were not computed in the IM-SRG(2) timings. 

For our second illustration of operator evolution, we 
consider the generalized COM Hamiltonian, Eq. In 
calculations of nuclei, the ground state expectation value 
of this quantity is useful to diagnose whether approxi¬ 
mate solutions of the Schrodinger equation are contami¬ 
nated by spurious COM motion. Since nuclei are self¬ 
bound objects governed by a translationally-invariant 
Hamiltonian, an exact solution of the Schrodinger equa¬ 
tion must factorize into the product of a wave function 
for the physically relevant intrinsic motion times a wave 
function for the COM coordinate, 

1^) = Win G Wcm 


(46) 
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FIG. 4. (Color online) Magnitude of the 0-body contributions of the term in Eq. 30 evaluated in the N02B approximation. 


The top row is for the electron gas at Wigner-Seitz radii of (a) = 0.5 and (b) Ts = 5.0. The bottom row is for ^°0, starting 

from the chiral NN potential of Entem and Machleidt [3H], softened by a free-space SRG evolution to (c) A = 2.0 fm“^ and (d) 
A = 3.0 fm~^. The electron gas calculations were done for = 14 electrons in a periodic box with M — 114 single particle 
orbitals. The ^®0 calculations were done in an Emax = 8 model space, with hui = 24 MeV for the underlying harmonic oscillator 
basis. 


As is well known, there are two strategies to rigor¬ 
ously guarantee this factorization; one can work in a 
translationally-invariant basis from the outset, or one can 
work in a so-called full NHut model space comprised of all 
A-particle harmonic oscillator Slater determinants with 
excitations up to and including Nhw. Neither choice is 
optimal since the former is limited to light nuclei due 
to the factorial scaling of the required antisymmetriza- 
tion, while the latter limits the choice of the single parti¬ 
cle orbitals to the harmonic oscillator basis and doesn’t 
carry over to methods that are capable of reaching heav¬ 
ier nuclei, such as coupled cluster theory and the IM-SRG 


where it is more natural to define the model space via an 
energy cutoff (e.g., 2n -|- 1 < Cmax) on the single particle 
states. In the case of calculations with an Cmax cutoff, 
there is no analytical guarantee that the COM and in¬ 
trinsic wave functions factorize. 

In Ref. [3^ , Hagen and collaborators gave an ingenious 
prescription to diagnose whether or not Eq. is satis¬ 
fied in such calculations. The basic idea is to assume 
that the factorized COM wave function is a Gaussian, 
and is therefore the ground state of Hcm(w) with eigen¬ 
value zero. Note that Co ^ ui in general, where oo is the 
frequency of the underlying oscillator basis. The pre- 
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FIG. 5. (Color online) Flowing IM-SRG(2) and Magnus(2) 
HEG correlation energy, Eq{s)—Ehf, for Wigner-Seitz radii of 
a) Vs = 5.0 and b) = 0.5. The solid black line corresponds 
to IM-SRG(2) results using an adaptive solver based on the 
Adams-Bashforth method, while the other lines correspond to 
Magnus(2) and IM-SRG(2) results using different Euler step 
sizes. The red circles denote the quasi-exact FCIQMC results 
of Ref. [35] . 


scription to find w involves solving a quadratic equation 


FIG. 6. (Color online) Flowing IM-SRG(2) and Magnus(2) 
ground state energy, Eo{s), for starting from the N^LO 
NN interaction of Entem and Machleidt |38| evolved by the 
free-space SRG to a) A = 2.7 fm“^ and A = 2.0 fm“^. The 
solid black line corresponds to IM-SRG(2) results using an 
adaptive solver based on the Adams-Bashforth method, while 
the other lines correspond to Magnus(2) and IM-SRG(2) re¬ 
sults using different Euler step sizes. The calculations were 
done in an Cmax = 8 model space, with hui = 24 MeV for the 
underlying harmonic oscillator basis. 


2/4 4 

hlb = hhj+ ± Y -(£'cm(w))^ + -hujEcmioj), 

(47) 

where 

E,M = (48) 

= lim (49) 

s—>-oo 

= . (50) 

s—^ 

Since there are two roots of Eq. we choose the one 
that gives a smaller value for i?cm(w). Applying this pre¬ 
scription to our calculations of ^®0, we obtain the results 
shown in Fig. In the top panel, we see that the expec¬ 
tation value of the COM Hamiltonian i7cm(t^) is approx¬ 
imately zero for w « 20 MeV, but varies parabolically 
and becomes rather large away from this point. This 
suggests that if Eq. [i^ is satisfied, the frequency of the 
factorized COM Gaussian should have w « 20 MeV. This 
is born out in the bottom panel, where the two roots of 


Eq. 47 are plotted as a function of hio. Picking the root 
that minimizes £'cm(0), we find that indeed w « 20 MeV 
over a wide range of w, and that ~ 0. Since 

the excitation energy for the first spurious COM mode 
is to « 20 MeV, while Ecm{^) ranges between 0.03-0.14 
MeV over the entire range of huj, we conclude that the 
factorization of COM/intrinsic motion is satisfied to an 
excellent approximation. 


V. SUMMARY AND OUTLOOK 


In this paper, we have shown how the Magnus expan¬ 
sion can be used to bypass computational limitations 
arising from the large memory demands of high-order 
ODE solvers typically used in IM-SRG calculations. The 
success of the Magnus expansion derives from the fact 
that by reformulating Eq. as a flow equation for the 
operator D(s), where U{s) = we are able to use a 

simple first-order Euler ODE solver without any loss of 
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FIG. 7. (Color online) (Color online) Electron gas momen¬ 
tum distributions calculated in the Magnus(2) approximation. 
The calculations were done for = 14 electrons in a periodic 
box with M — 778 single particle orbitals. 



FIG. 8. (Color online) Timing for Magnus(2) and IM-SRG(2) 
HEG calculations as the single particle basis is increased. The 
two bottom curves are for Vs = .5 and the top for = 5. For 
the Magnus(2) timing, this includes the calculations of the 
momentum distributions. 


accuracy, resulting in substantial memory savings and a 
modest reduction in CPU time. In conventional formu¬ 
lations of the SRG, time step errors accumulate directly 
in the evolved H{s), necessitating the use of a high-order 
solver to preserve an acceptable level of accuracy. In 
the Magnus expansion, even though sizable time step er¬ 
rors accumulate in n(s) with a first-order Euler method, 
upon exponentiation the transformation is still unitary, 
and the transformed H{s) = W {s)HU{s) is unitarily 



hco [MeV] 


FIG. 9. (Color online) (Color online) Center of mass diagnos¬ 
tics for Magnus(2) calculations of starting from the N^LO 
NN interaction of Entem and Machleidt |38| evolved by the 
free-space SRG to A = 2.0 fm“^. See the text for details. The 
calculations were done in an Cmax = 9 model space. 


equivalent to the initial Hamiltonian modulo any trunca¬ 
tions made (e.g., the N02B approximation or numerical 
truncation of the infinite series) in evaluating the BCH 
formula. 


After introducing the basic formalism of the Mag¬ 
nus expansion and illustrating its numerical virtues 
for a schematic matrix model, we turned to realistic 
many-body calculations of the homogeneous electron 
gas (HEG) and To make progress, we introduced 
the Magnus(2) approximation in which all operators 
(i/(s), B(s), ? 7 (s)) and commutators are truncated at the 
N02B approximation, and the non-terminating Magnus 
derivative, Eq. and BCH formula, Eq. are trun¬ 
cated numerically. In all cases studied, our calculations 
converge to a final answer that is independent of step 
size and agrees well with IM-SRG(2) calculations using 
a high-order solver. Moreover, our final results are in¬ 
dependent of the precise convergence criteria for trun¬ 
cating Eqs. 28 and ^ provided that Cderiv ^ 10“^ and 


ebch ^ 10 ^MeV (Hartree) for the (HEG) calcula¬ 
tions. 


The evaluation of observables besides the Hamiltonian 
poses considerable challenges in the IM-SRG since the 
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evolution of each additional observable roughly doubles 
the dimensionality of the flow equations. In the Mag¬ 
nus expansion formulation, the evolution of additional 
operators is trivial since one solves flow equations for the 
unitary transformation and then constructs the evolved 
observables by application of the BCH formula. The di¬ 
mensionality of the flow equations is fixed, regardless of 
how many additional operators are being evolved. Proof- 
of-principle operator evolutions were carried out using 
the Magnus expansion for the momentum distribution 
in the HEG, and the generalized COM Hamiltonian in 
This opens the door for the ab-initio calculation 
of a variety of properties (radii, transition matrix ele¬ 


ments, response functions, etc.) in addition to energies 
for medium-mass nuclei. 
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